#include "Header.h"
#include <cstdio>
#include <unistd.h>
int main(){
	/*  Test the basic benchmark function */
	Benchmarks* fp=NULL;
	int funToRun[] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};//function set
	int funNum = 20;//the total number of functions to be optimized in CEC'2010 benchmark set
    unsigned run_index;
    int i, j, t, k,TS1,TS2, e1, e2,NP;
    int index,times,mark;
	int gbest=0;//global best index
    int function_index;//functional index
    int fitness_counter;//Current fitness evaluation times
    double final_val;// the final global best solution
    double MAX;
    double MIN;
    double Beta;//learning rate
    double Wmax;//maximum weight
    double Wmin;//minimum weight
    double **population = new double*[NPmax];//store the positions of the population
	double **speed = new double*[NPmax];//store the velocity of the population
    double **position = new double*[NPmax];//store the positions of the population
	double **velocity = new double*[NPmax];//store the velocity of the population
    double r,R1, R2, R3;//random number
    double *results = new double[NPmax];// the fitness results of the whole population
    double *value = new double[NPmax];
    double *w = new double[NPmax];//weight
    double *p = new double[NPmax];//probability
    int *selected_particles = new int[NPmax];
    struct fitnessIndexPairs
    {
        double fitnessValue;
        int particleIndex;
    };
    fitnessIndexPairs* particles = new fitnessIndexPairs[NPmax];
    for (i = 0; i < NPmax; ++i)
	{
		population[i] = new double[dim];
		speed[i] = new double[dim];
		position[i] = new double[dim];
		velocity[i] = new double[dim];
	}
    cout << "Start to optimize the function!" << endl;
    for (function_index = 0; function_index < funNum; function_index++)
    {
        cout << "Function " << funToRun[function_index] << " Begined!" << endl;
        fp = generateFuncObj(funToRun[function_index]);
        MIN = fp->getMinX();
        MAX = fp->getMaxX();
        boost::mt19937 generator(time(0)*rand());
        boost::uniform_real<> uniform_real_generator2(MIN, MAX);
        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random(generator, uniform_real_generator2);
        boost::uniform_real<> uniform_real_generate_r(0, 1);
        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_num_r(generator, uniform_real_generate_r);//to generate a random number within [0,1]
        // Gaussian distribution for generating the learning rate Beta
        boost::normal_distribution<double> normal_distribution1(0.4, 0.01);
        boost::variate_generator<boost::mt19937&, boost::normal_distribution<double>> normal_distribution2(generator, normal_distribution1);
        boost::uniform_int<> uniform_int_generator( TSmin, TSmax);
        boost::variate_generator< boost::mt19937&, boost::uniform_int<> > random_select_TS( generator, uniform_int_generator );
        char fun[20];
        snprintf(fun, 20, "%d", funToRun[function_index]);
        string filename_fitness = "./results/Fitness_result_for_"+ string(fun) + ".txt";
        string filename_fitness1 = "./results/fv/Fitness_result_for_"+ string(fun) + ".txt";
        ofstream out_fitness(filename_fitness.c_str());
        ofstream out_fitness1(filename_fitness1.c_str());
        if (!out_fitness)
        {
            cerr << "Can not open the file " << filename_fitness << endl;
            exit(1);
        }
        for (run_index = 1; run_index <= timesOfRun; run_index++)
        {
            cout << "Running the " << run_index << "th time!!" << endl;
            fitness_counter = 0;
			// to initialize the population
			for (i = 0; i < NPmax; i++)
			{
				for (j = 0; j < dim; j++)
				{
					population[i][j] = random();
					position[i][j]=population[i][j];
					speed[i][j] = 0;
					velocity[i][j]=speed[i][j];
				}
				results[i] = fp->compute(population[i]);
				value[i]=results[i];
				particles[i].fitnessValue = results[i];
                particles[i].particleIndex = i;
			}
			final_val = results[0];
            for (i = 1; i < NPmax; i++)
            {
                if (final_val > results[i])
                {
                    final_val = results[i];
                    gbest=i;
                }
            }
            fitness_counter += NPmax;
            times=0;
            NP=NPmax;
            //main iterations
            while (fitness_counter < MAX_FV)
            {
                // sort particle according to its fitness
                std::sort(particles, particles + NP, [](const fitnessIndexPairs& a, const fitnessIndexPairs& b) {
                    return a.fitnessValue < b.fitnessValue; // from best to worst
                });
                Wmax=(1/ (2 * NP*sqrt(2 * PI)));
                Wmin=(1/ (2 * NP*sqrt(2 * PI))) * exp(-(pow((NP-1), 2)/ (8 *PI* pow(NP, 2))));
                for (i = 0; i < NP; ++i)
                {
                    w[i]=(1/ (2 * NP*sqrt(2 * PI))) * exp(-(pow(i, 2) / (8 *PI* pow(NP, 2))));
                    p[i]=(Wmax - w[i])/(Wmax - Wmin);
                }
                //counter=0;
                for (i = NP-1; i > TSmax; i--)
                {
                    r = random_real_num_r();
                    if(r<p[i])
                    {
                        TS1 = random_select_TS();
                        TS2 = random_select_TS();
                        for (t = 0; t < i; t++)
                        {
                            selected_particles[t] = t;
                        }
                        std::shuffle(selected_particles, selected_particles + i, generator);
                        //find the best particle Xe1 in TS1
                        e1=selected_particles[0];
                        for(t = 1; t < TS1; t++)
                        {
                            if(selected_particles[t]<e1)
                            {
                                e1=selected_particles[t];
                            }
                        }
                        //find another best particle Xe2 in TS2
                        do
                        {
                            std::shuffle(selected_particles, selected_particles + i, generator);
                            e2=selected_particles[0];
                            for(t = 1; t < TS2; t++)
                            {
                                if(selected_particles[t]<e2)
                                {
                                    e2=selected_particles[t];
                                }
                            }
                        }while(particles[e2].particleIndex==particles[e1].particleIndex);
                        // Check if f(e1) < f(e2)
                        if (particles[e2].fitnessValue < particles[e1].fitnessValue)
                        {
                            // Swap e1 and e2
                            swap( e1, e2 );
                        }
                        do {
                                Beta = normal_distribution2();//generate the Beta vale
                            } while (Beta <= 0.0 || Beta >= 1.0); //Bata is set within[0,1]
                        //update particle
                        for (j = 0; j < dim; ++j)
                        {
                            //randomly generate two random numbers
                            R1 = random_real_num_r();
                            R2 = random_real_num_r();
                            R3 = random_real_num_r();
                            //update the velocity
                            speed[particles[i].particleIndex][j] = R1 * speed[particles[i].particleIndex][j] + R2 * (population[particles[e1].particleIndex][j] - population[particles[i].particleIndex][j]) + Beta * R3 * (population[particles[e2].particleIndex][j] - population[particles[i].particleIndex][j]);
                            velocity[particles[i].particleIndex][j]=speed[particles[i].particleIndex][j];
                            //update the position
                            population[particles[i].particleIndex][j] = population[particles[i].particleIndex][j] + speed[particles[i].particleIndex][j];
                            //check whether the position is out of the range
                            if (population[particles[i].particleIndex][j] < MIN)
                            {
                                population[particles[i].particleIndex][j] = MIN;
                            }
                            if (population[particles[i].particleIndex][j] > MAX)
                            {
                                population[particles[i].particleIndex][j] = MAX;
                            }
                            position[particles[i].particleIndex][j]=population[particles[i].particleIndex][j];
                        }
                        particles[i].fitnessValue = fp->compute(population[particles[i].particleIndex]);
                        results[particles[i].particleIndex] = particles[i].fitnessValue;
                        value[particles[i].particleIndex]=results[particles[i].particleIndex];
                        if (final_val > results[particles[i].particleIndex])
                        {
                            final_val = results[particles[i].particleIndex];
                            gbest=particles[i].particleIndex;
                        }
                        fitness_counter++;
                    }
                }
                k=NP;
                //calculate the current NP
                NP=ceil(NPmax+(NPmin-NPmax)*(double)(fitness_counter)/(double)(MAX_FV));
                if(NP<k)//if NP becomes smaller
                {
                    //delete superfluous particles except the best one
                    index = 0;
                    for (t = 0; t < k; t++)
                    {
                        if (t != gbest)
                        {
                            selected_particles[index] = t;
                            index++;
                        }
                    }
                    std::shuffle(selected_particles, selected_particles + (k - 1), generator);
                    for (j = 0; j < dim; ++j)
                    {
                        population[0][j]=position[gbest][j];
                        speed[0][j]=velocity[gbest][j];
                    }
                    results[0]=value[gbest];
                    particles[0].fitnessValue = results[0];
                    particles[0].particleIndex = 0;
                    mark=0;
                    for (t = 1; t < NP; t++)
                    {
                        for (j = 0; j < dim; ++j)
                        {
                            population[t][j]=position[selected_particles[mark]][j];
                            speed[t][j]=velocity[selected_particles[mark]][j];
                        }
                        results[t]=value[selected_particles[mark]];
                        particles[t].fitnessValue = results[t];
                        particles[t].particleIndex = t;
                        mark++;
                    }
                }
                if (fitness_counter >= 100000*times)
                {
                    out_fitness1 << final_val << "\t";
                    times++;
                }
            }//end of while
            out_fitness1 << endl;
            out_fitness << final_val << endl;  // Move to the next line for the next independent experiment
            //printf("F %d value = %1.20E\n", fp->getID(), final_val);
        }//end of for runs
        cout << "Function " << funToRun[function_index] << " Finished!" << endl;
        delete fp;
        out_fitness.close();
        out_fitness1.close();
    }
    //release the resources
    for (t = 0; t < NPmax; ++t)
    {
        delete[]population[t];
        delete[]speed[t];
        delete[]position[t];
	delete[]velocity[t];
    }
    delete[]population;
    delete[]speed;
    delete[]position;
    delete[]velocity;
    delete[]results;
    delete[]value;
    delete[]w;
    delete[]p;
    delete[]selected_particles;
    delete[]particles;
    return 0;
}
// create new object of class with default setting
Benchmarks* generateFuncObj(int funcID){
	Benchmarks *fp;
	// run each of specified function in "configure.ini"
	if (funcID==1){
		fp = new F1();
	}else if (funcID==2){
		fp = new F2();
	}else if (funcID==3){
		fp = new F3();
	}else if (funcID==4){
		fp = new F4();
	}else if (funcID==5){
		fp = new F5();
	}else if (funcID==6){
		fp = new F6();
	}else if (funcID==7){
		fp = new F7();
	}else if (funcID==8){
		fp = new F8();
	}else if (funcID==9){
		fp = new F9();
	}else if (funcID==10){
		fp = new F10();
	}else if (funcID==11){
		fp = new F11();
	}else if (funcID==12){
		fp = new F12();
	}else if (funcID==13){
		fp = new F13();
	}else if (funcID==14){
		fp = new F14();
	}else if (funcID==15){
		fp = new F15();
	}else if (funcID==16){
		fp = new F16();
	}else if (funcID==17){
		fp = new F17();
	}else if (funcID==18){
		fp = new F18();
	}else if (funcID==19){
		fp = new F19();
	}else if (funcID==20){
		fp = new F20();
	}else{
		cerr<<"Fail to locate Specified Function Index"<<endl;
		exit(-1);
	}
	return fp;
}
